Fast track to GRASS with R: the rgrass package

In the same way there are Python packages that allow us to interact with GRASS GIS tools and data, there is an R package, rgrass, that serves as the interface between GRASS and R. The rgrass package is developed and maintained by @rgrass and can be found at: https://github.com/rsbivand/rgrass/. In this fast track tutorial, we’ll learn how to use GRASS GIS from R.

Main funtions

The main functions within rgrass are the following:

initGRASS(): it starts a GRASS GIS session from R.
execGRASS(): it executes GRASS GIS commands from R.
gmeta(): it prints GRASS GIS session metadata: database, project, mapset, computational region settings and CRS details.
read_VECT() and read_RAST(): they read vector and raster maps from a GRASS project into R.
write_VECT() and write_RAST(): they write vector and raster objects from R into a GRASS GIS project.

Note

For further details on rgrass functionality, use examples and data format coercion, see: https://rsbivand.github.io/rgrass/.

Basic Usage: Choose your own adventure

If you already use GRASS as your geospatial data processing engine, you most likely have your spatial data within GRASS projects. You might need however to do some statistical analysis, some modelling and prediction or create publication ready visualizations in R. In such cases, similar to how it works for the GRASS Python API within Jupyter notebooks (), you can start a GRASS session from R or RStudio.

If you are not a regular GRASS user, you probably don’t care about GRASS projects, and you just want to use GRASS functionality because you know it rocks! Well, rgrass gets you covered in such cases, too!

Let’s see the general basic steps and then we’ll dive into the details:

Open R (or RStudio)
Load rgrass library with library(rgrass)
Start a GRASS GIS session with initGRASS() in an existent or temporary project
Use GRASS GIS modules through execGRASS()
Use read_VECT(), read_RAST(), write_VECT() and write_RAST() to read data from and write data into GRASS database

A. You have GRASS projects

Let’s see an example for the case when we do all our geospatial data processing within GRASS and hence have all the spatial data organized within projects and mapsets but we need to run some statistical analysis, modelling and prediction or visualization in R.

We start R or Rstudio and load the rgrass library. It will tell us that GRASS is not running, but we know that already… and that’s about to change in a moment.

library(rgrass)

GRASS GIS interface loaded with GRASS version: (GRASS not running)

We’ll next start GRASS from within R or RStudio using the initGRASS() function. Since we want to start GRASS in a specific project and mapset, we need to specify them. Optionally, we specify which GRASS binary to use. This might be useful in the case we have many GRASS versions in our system. If not provided, initGRASS() will attempt to find it in default locations depending on your operative system.

We can now list and read our GRASS raster and vector maps into R and do our statistical analysis, modelling and/or visualizations using other R packages. Obviously, we can also use rgrass as a mere interface to analyze data within GRASS, i.e., no reading from or writing to GRASS needed, we can just use GRASS from R, equivalent to what we do with the GRASS Python API (). Here, we’ll demonstrate the use of all the main rgrass functions mentioned above.

Let’s then list our GRASS raster and vector maps:

# list GRASS raster maps (equivalent to `g.list type=raster`)
execGRASS("g.list", parameters = list(type="raster"))

basins
developed
elevation
elevation_shade
geology
lakes
landuse
result_from_R
soils

# list GRASS vector maps (equivalent to `g.list type=vector`)
execGRASS("g.list", parameters = list(type="vector"))

boundary_region
boundary_state
census
elev_points
firestations
geology
geonames
hospitals
points_of_interest
railroads
roadsmajor
schools
streams
streets
zipcodes

The resulting map lists could be saved in an R object that we can subset later in case we want to import several but not all raster maps, for example. Let’s see how to do that.

# save map list in an object
rast_list <- execGRASS("g.list", parameters = list(type="raster"))

basins
developed
elevation
elevation_shade
geology
lakes
landuse
result_from_R
soils

rast_list

[1] 0
attr(,"resOut")
[1] "basins"          "developed"       "elevation"       "elevation_shade"
[5] "geology"         "lakes"           "landuse"         "result_from_R"  
[9] "soils"          
attr(,"resErr")
character(0)

# retrieve only the maps list and overwrite rast_list
rast_list <- attributes(rast_list)$resOut

# import elevation and landuse
to_import <- rast_list[c(3,7)]

maplist <- list()
for (i in to_import) {
  maplist[i] <- read_RAST(i)
}

Checking GDAL data type and nodata value...
   2%   5%   8%  11%  14%  17%  20%  23%  26%  29%  32%  35%  38%  41%  44%  47%  50%  53%  56%  59%  62%  65%  68%  71%  74%  77%  80%  83%  86%  89%  92%  95%  98% 100%
Using GDAL data type <Float32>
Exporting raster data to RRASTER format...
   2%   5%   8%  11%  14%  17%  20%  23%  26%  29%  32%  35%  38%  41%  44%  47%  50%  53%  56%  59%  62%  65%  68%  71%  74%  77%  80%  83%  86%  89%  92%  95%  98% 100%
r.out.gdal complete. File </tmp/Rtmp6pJLgV/file218d3396e8f3a.grd> created.

Warning in `[<-`(`*tmp*`, i, value = read_RAST(i)): implicit list embedding of
S4 objects is deprecated

Checking GDAL data type and nodata value...
   2%   5%   8%  11%  14%  17%  20%  23%  26%  29%  32%  35%  38%  41%  44%  47%  50%  53%  56%  59%  62%  65%  68%  71%  74%  77%  80%  83%  86%  89%  92%  95%  98% 100%
Using GDAL data type <UInt32>
Exporting raster data to RRASTER format...
   2%   5%   8%  11%  14%  17%  20%  23%  26%  29%  32%  35%  38%  41%  44%  47%  50%  53%  56%  59%  62%  65%  68%  71%  74%  77%  80%  83%  86%  89%  92%  95%  98% 100%
r.out.gdal complete. File </tmp/Rtmp6pJLgV/file218d34740771e.grd> created.

Warning in `[<-`(`*tmp*`, i, value = read_RAST(i)): implicit list embedding of
S4 objects is deprecated

maplist

$elevation
class       : SpatRaster 
dimensions  : 1350, 1500, 1  (nrow, ncol, nlyr)
resolution  : 10, 10  (x, y)
extent      : 630000, 645000, 215000, 228500  (xmin, xmax, ymin, ymax)
coord. ref. : NAD83(HARN) / North Carolina (EPSG:3358) 
source      : file218d3396e8f3a.grd 
name        : file218d3396e8f3a 
min value   :          55.57879 
max value   :         156.32986 

$landuse
class       : SpatRaster 
dimensions  : 1350, 1500, 1  (nrow, ncol, nlyr)
resolution  : 10, 10  (x, y)
extent      : 630000, 645000, 215000, 228500  (xmin, xmax, ymin, ymax)
coord. ref. : NAD83(HARN) / North Carolina (EPSG:3358) 
source      : file218d34740771e.grd 
name        : file218d34740771e 
min value   :                 0 
max value   :                 7

GRASS raster and vector maps are translated into terra’s packege SpatRaster and SpatVector objects, respectively. These objects can then, within R, be easily coerced to other types of spatial objects such as simple features (sf), stars, etc.

Remember that raster objects will always be exported from GRASS following the computational region settings. So, bear that in mind when reading into R which will hold them in memory. Vectors however will be exported in their full extent.

Note

See terra vignettes with further explanations and examples: https://rspatial.github.io/terra/.

Let’s load the terra library to quickly display our recently imported raster maps:

library(terra)

terra 1.7.71

plot(maplist$elevation)

Optionally, we could stack our two SpatRaster objects together and plot them together:

rstack <- rast(maplist)
plot(rstack)

Let’s create a boxplot of elevation per land class.

boxplot(rstack$elevation, rstack$landuse, maxcell=50000)

Warning: [boxplot] taking a regular sample of 50000 cells

Let’s import a vector map, too, and explore its attributes.

census <- read_VECT("census")

Warning in execGRASS("v.out.ogr", flags = flags, input = vname, type = type, : The command:
v.out.ogr --overwrite input=census type=area layer=1 output=/tmp/Rtmp6pJLgV/file218d31e29fc36.gpkg output_layer=census format=GPKG
produced at least one warning during execution:
Exporting 2547 areas (may take some time)...
   5%  11%  17%  23%  29%  35%  41%  47%  53%  59%  65%  71%  77%  83%  89%  95% 100%
WARNING: 10 features without category were skipped. Features without
         category are written only when -c flag is given.
v.out.ogr complete. 2537 features (Polygon type) written to <census> (GPKG
format).

Exporting 2547 areas (may take some time)...
   5%  11%  17%  23%  29%  35%  41%  47%  53%  59%  65%  71%  77%  83%  89%  95% 100%
WARNING: 10 features without category were skipped. Features without
         category are written only when -c flag is given.
v.out.ogr complete. 2537 features (Polygon type) written to <census> (GPKG
format).

head(census)

  cat OBJECTID      AREA PERIMETER BLOCK_ BLOCK_ID        BLOCKNUM TOTAL_POP
1  86    82802  9981.633   444.547  82803    82802 371830535104051        13
2  88    82805   304.076    91.371  82806    82805 371830535104046         0
3  91    82811  2400.879   192.304  82812    82811 371830535104034        20
4 140    82955  1239.619   147.011  82956    82955 371830535104039         0
5 173    83073 10358.715   429.738  83074    83073 371830535014007        80
6 228    83352 15459.671   458.961  83353    83352 371830535015007        18
  POP_1RACE WHITE_ONLY BLACK_ONLY AMIND_ONLY ASIAN_ONLY HWPAC_ONLY OTHER_ONLY
1        13         13          0          0          0          0          0
2         0          0          0          0          0          0          0
3        20          3          0          0         17          0          0
4         0          0          0          0          0          0          0
5        76         38          9          0          0          0         29
6        18         17          1          0          0          0          0
  POP_2RACES HISPANIC MALE FEMALE P_UNDER_5 P5_TO_9 P10_TO_14 P15_TO_19
1          0        0    6      7         1       2         1         1
2          0        0    0      0         0       0         0         0
3          0        0   10     10         1       1         1         1
4          0        0    0      0         0       0         0         0
5          4       54   37     43        20       4         5         7
6          0        3   10      8         1       0         1         0
  P20_TO_24 P25_TO_34 P35_TO_44 P45_TO_54 P55_TO_59 P60_TO_64 P65_TO_69
1         1         3         1         1         1         0         1
2         0         0         0         0         0         0         0
3         2         8         2         3         0         1         0
4         0         0         0         0         0         0         0
5        14        17         7         5         0         0         0
6         0         5         0         5         1         2         2
  P70_TO_74 P75_TO_84 P75_OLDER P85_OLDER MEDIAN_AGE P18_OLDER P21_OLDER
1         0         0         0         0       25.5         9         8
2         0         0         0         0        0.0         0         0
3         0         0         0         0       30.0        16        16
4         0         0         0         0        0.0         0         0
5         0         1         1         0       21.8        48        44
6         0         1         1         0       52.5        16        16
  P65_OLDER M75_OLDER F75_OLDER HOUSEHOLDS HH_SIZE SINGLE_HH FAMILY_HH
1         1         0         0          3    4.33         0         2
2         0         0         0          0    0.00         0         0
3         0         0         0          7    2.86         0         6
4         0         0         0          0    0.00         0         0
5         1         0         1         17    4.71         1        16
6         3         0         1          9    2.00         3         6
  MAR_COUPLE W_OWN_CHLD M_SNGL_PAR F_SNGL_PAR NONFAM_HH HH_W_CHILD HH_W_ELDER
1          2          2          0          0         1          3          1
2          0          0          0          0         0          0          0
3          6          3          0          0         1          3          0
4          0          0          0          0         0          0          0
5          5          3          3          3         0         10          1
6          6          2          0          0         0          2          2
  SNGL_ELDER POP_IN_HH FAM_HHLDER SPOUSE CHILD GROUP_QTRS GQ_INST GQ_NONINST
1          0        13          2      2     4          0       0          0
2          0         0          0      0     0          0       0          0
3          0        20          6      6     6          0       0          0
4          0         0          0      0     0          0       0          0
5          0        80         16      5    22          0       0          0
6          1        18          6      6     3          0       0          0
  FAM_SIZE HOUSING_U OCCUPIED VACANT OWNER_U RENTER_U VACANT_SEA HH_SIZE_OW
1     4.00         3        3      0       0        3          0       0.00
2     0.00         0        0      0       0        0          0       0.00
3     3.00         9        7      2       7        0          0       2.86
4     0.00         0        0      0       0        0          0       0.00
5     3.81        17       17      0       9        8          0       5.00
6     2.50         9        9      0       8        1          0       1.88
  HH_SIZE_RE
1          4
2          0
3          0
4          0
5          4
6          3

summary(census)

      cat          OBJECTID           AREA           PERIMETER       
 Min.   :   1   Min.   : 82370   Min.   :    155   Min.   :   60.03  
 1st Qu.: 635   1st Qu.:136909   1st Qu.:  11924   1st Qu.:  489.03  
 Median :1269   Median :140048   Median :  25677   Median :  728.51  
 Mean   :1269   Mean   :130677   Mean   :  91514   Mean   : 1134.20  
 3rd Qu.:1903   3rd Qu.:143047   3rd Qu.:  65508   3rd Qu.: 1268.41  
 Max.   :2537   Max.   :149473   Max.   :5392753   Max.   :14045.55  
     BLOCK_          BLOCK_ID        BLOCKNUM           TOTAL_POP      
 Min.   : 82371   Min.   : 82370   Length:2537        Min.   :   0.00  
 1st Qu.:136910   1st Qu.:136909   Class :character   1st Qu.:   0.00  
 Median :140049   Median :140048   Mode  :character   Median :  26.00  
 Mean   :130678   Mean   :130677                      Mean   :  63.03  
 3rd Qu.:143048   3rd Qu.:143047                      3rd Qu.:  62.00  
 Max.   :149474   Max.   :149473                      Max.   :6173.00  
   POP_1RACE         WHITE_ONLY        BLACK_ONLY        AMIND_ONLY     
 Min.   :   0.00   Min.   :   0.00   Min.   :   0.00   Min.   : 0.0000  
 1st Qu.:   0.00   1st Qu.:   0.00   1st Qu.:   0.00   1st Qu.: 0.0000  
 Median :  26.00   Median :   9.00   Median :   0.00   Median : 0.0000  
 Mean   :  61.85   Mean   :  40.02   Mean   :  17.18   Mean   : 0.2598  
 3rd Qu.:  61.00   3rd Qu.:  43.00   3rd Qu.:  12.00   3rd Qu.: 0.0000  
 Max.   :6054.00   Max.   :4645.00   Max.   :1138.00   Max.   :55.0000  
   ASIAN_ONLY        HWPAC_ONLY        OTHER_ONLY        POP_2RACES     
 Min.   :  0.000   Min.   : 0.0000   Min.   :  0.000   Min.   :  0.000  
 1st Qu.:  0.000   1st Qu.: 0.0000   1st Qu.:  0.000   1st Qu.:  0.000  
 Median :  0.000   Median : 0.0000   Median :  0.000   Median :  0.000  
 Mean   :  2.402   Mean   : 0.0335   Mean   :  1.953   Mean   :  1.179  
 3rd Qu.:  0.000   3rd Qu.: 0.0000   3rd Qu.:  0.000   3rd Qu.:  0.000  
 Max.   :401.000   Max.   :16.0000   Max.   :243.000   Max.   :119.000  
    HISPANIC            MALE             FEMALE          P_UNDER_5      
 Min.   :  0.000   Min.   :   0.00   Min.   :   0.00   Min.   :  0.000  
 1st Qu.:  0.000   1st Qu.:   0.00   1st Qu.:   0.00   1st Qu.:  0.000  
 Median :  0.000   Median :  12.00   Median :  13.00   Median :  1.000  
 Mean   :  4.036   Mean   :  31.88   Mean   :  31.15   Mean   :  3.534  
 3rd Qu.:  1.000   3rd Qu.:  31.00   3rd Qu.:  31.00   3rd Qu.:  3.000  
 Max.   :704.000   Max.   :3714.00   Max.   :2459.00   Max.   :143.000  
    P5_TO_9          P10_TO_14        P15_TO_19          P20_TO_24      
 Min.   :  0.000   Min.   :  0.00   Min.   :   0.000   Min.   :   0.00  
 1st Qu.:  0.000   1st Qu.:  0.00   1st Qu.:   0.000   1st Qu.:   0.00  
 Median :  1.000   Median :  1.00   Median :   1.000   Median :   1.00  
 Mean   :  3.595   Mean   :  3.47   Mean   :   5.594   Mean   :   8.99  
 3rd Qu.:  3.000   3rd Qu.:  3.00   3rd Qu.:   3.000   3rd Qu.:   4.00  
 Max.   :127.000   Max.   :112.00   Max.   :3419.000   Max.   :2499.00  
   P25_TO_34        P35_TO_44        P45_TO_54         P55_TO_59     
 Min.   :  0.00   Min.   :  0.00   Min.   :  0.000   Min.   : 0.000  
 1st Qu.:  0.00   1st Qu.:  0.00   1st Qu.:  0.000   1st Qu.: 0.000  
 Median :  3.00   Median :  4.00   Median :  3.000   Median : 1.000  
 Mean   : 11.49   Mean   :  9.61   Mean   :  7.443   Mean   : 2.349  
 3rd Qu.: 11.00   3rd Qu.: 10.00   3rd Qu.:  8.000   3rd Qu.: 3.000  
 Max.   :467.00   Max.   :404.00   Max.   :180.000   Max.   :68.000  
   P60_TO_64        P65_TO_69        P70_TO_74        P75_TO_84     
 Min.   : 0.000   Min.   : 0.000   Min.   : 0.000   Min.   : 0.000  
 1st Qu.: 0.000   1st Qu.: 0.000   1st Qu.: 0.000   1st Qu.: 0.000  
 Median : 0.000   Median : 0.000   Median : 0.000   Median : 0.000  
 Mean   : 1.768   Mean   : 1.416   Mean   : 1.331   Mean   : 1.843  
 3rd Qu.: 2.000   3rd Qu.: 2.000   3rd Qu.: 2.000   3rd Qu.: 2.000  
 Max.   :40.000   Max.   :55.000   Max.   :37.000   Max.   :81.000  
   P75_OLDER         P85_OLDER          MEDIAN_AGE      P18_OLDER      
 Min.   :  0.000   Min.   :  0.0000   Min.   : 0.00   Min.   :   0.00  
 1st Qu.:  0.000   1st Qu.:  0.0000   1st Qu.: 0.00   1st Qu.:   0.00  
 Median :  0.000   Median :  0.0000   Median :30.80   Median :  21.00  
 Mean   :  2.442   Mean   :  0.5987   Mean   :25.93   Mean   :  50.49  
 3rd Qu.:  3.000   3rd Qu.:  0.0000   3rd Qu.:40.10   3rd Qu.:  51.00  
 Max.   :220.000   Max.   :139.0000   Max.   :85.60   Max.   :6147.00  
   P21_OLDER         P65_OLDER         M75_OLDER         F75_OLDER      
 Min.   :   0.00   Min.   :  0.000   Min.   : 0.0000   Min.   :  0.000  
 1st Qu.:   0.00   1st Qu.:  0.000   1st Qu.: 0.0000   1st Qu.:  0.000  
 Median :  20.00   Median :  2.000   Median : 0.0000   Median :  0.000  
 Mean   :  44.69   Mean   :  5.189   Mean   : 0.8065   Mean   :  1.635  
 3rd Qu.:  49.00   3rd Qu.:  6.000   3rd Qu.: 1.0000   3rd Qu.:  2.000  
 Max.   :1429.00   Max.   :230.000   Max.   :43.0000   Max.   :177.000  
   HOUSEHOLDS        HH_SIZE        SINGLE_HH         FAMILY_HH     
 Min.   :  0.00   Min.   :0.000   Min.   :  0.000   Min.   :  0.00  
 1st Qu.:  0.00   1st Qu.:0.000   1st Qu.:  0.000   1st Qu.:  0.00  
 Median : 11.00   Median :2.000   Median :  2.000   Median :  5.00  
 Mean   : 24.14   Mean   :1.617   Mean   :  7.695   Mean   : 13.11  
 3rd Qu.: 27.00   3rd Qu.:2.520   3rd Qu.:  8.000   3rd Qu.: 14.00  
 Max.   :663.00   Max.   :7.000   Max.   :281.000   Max.   :360.00  
   MAR_COUPLE        W_OWN_CHLD        M_SNGL_PAR        F_SNGL_PAR   
 Min.   :  0.000   Min.   :  0.000   Min.   : 0.0000   Min.   : 0.00  
 1st Qu.:  0.000   1st Qu.:  0.000   1st Qu.: 0.0000   1st Qu.: 0.00  
 Median :  3.000   Median :  1.000   Median : 0.0000   Median : 0.00  
 Mean   :  9.358   Mean   :  4.158   Mean   : 0.3366   Mean   : 1.76  
 3rd Qu.: 10.000   3rd Qu.:  4.000   3rd Qu.: 0.0000   3rd Qu.: 1.00  
 Max.   :225.000   Max.   :150.000   Max.   :13.0000   Max.   :79.00  
   NONFAM_HH         HH_W_CHILD        HH_W_ELDER        SNGL_ELDER     
 Min.   :  0.000   Min.   :  0.000   Min.   :  0.000   Min.   :  0.000  
 1st Qu.:  0.000   1st Qu.:  0.000   1st Qu.:  0.000   1st Qu.:  0.000  
 Median :  0.000   Median :  2.000   Median :  1.000   Median :  0.000  
 Mean   :  3.342   Mean   :  6.864   Mean   :  3.696   Mean   :  1.566  
 3rd Qu.:  3.000   3rd Qu.:  7.000   3rd Qu.:  5.000   3rd Qu.:  2.000  
 Max.   :178.000   Max.   :236.000   Max.   :158.000   Max.   :146.000  
   POP_IN_HH         FAM_HHLDER         SPOUSE            CHILD       
 Min.   :   0.00   Min.   :  0.00   Min.   :  0.000   Min.   :  0.00  
 1st Qu.:   0.00   1st Qu.:  0.00   1st Qu.:  0.000   1st Qu.:  0.00  
 Median :  25.00   Median :  5.00   Median :  3.000   Median :  4.00  
 Mean   :  56.56   Mean   : 13.11   Mean   :  9.358   Mean   : 13.77  
 3rd Qu.:  61.00   3rd Qu.: 14.00   3rd Qu.: 10.000   3rd Qu.: 14.00  
 Max.   :1676.00   Max.   :360.00   Max.   :225.000   Max.   :461.00  
   GROUP_QTRS          GQ_INST           GQ_NONINST          FAM_SIZE    
 Min.   :   0.000   Min.   :   0.000   Min.   :   0.000   Min.   :0.000  
 1st Qu.:   0.000   1st Qu.:   0.000   1st Qu.:   0.000   1st Qu.:0.000  
 Median :   0.000   Median :   0.000   Median :   0.000   Median :2.500  
 Mean   :   6.469   Mean   :   1.734   Mean   :   4.734   Mean   :1.913  
 3rd Qu.:   0.000   3rd Qu.:   0.000   3rd Qu.:   0.000   3rd Qu.:3.000  
 Max.   :5952.000   Max.   :1222.000   Max.   :5952.000   Max.   :7.500  
   HOUSING_U        OCCUPIED          VACANT           OWNER_U      
 Min.   :  0.0   Min.   :  0.00   Min.   :  0.000   Min.   :  0.00  
 1st Qu.:  0.0   1st Qu.:  0.00   1st Qu.:  0.000   1st Qu.:  0.00  
 Median : 12.0   Median : 11.00   Median :  0.000   Median :  4.00  
 Mean   : 25.8   Mean   : 24.14   Mean   :  1.653   Mean   : 12.16  
 3rd Qu.: 28.0   3rd Qu.: 27.00   3rd Qu.:  2.000   3rd Qu.: 14.00  
 Max.   :729.0   Max.   :663.00   Max.   :115.000   Max.   :280.00  
    RENTER_U        VACANT_SEA        HH_SIZE_OW       HH_SIZE_RE   
 Min.   :  0.00   Min.   : 0.0000   Min.   : 0.000   Min.   :0.000  
 1st Qu.:  0.00   1st Qu.: 0.0000   1st Qu.: 0.000   1st Qu.:0.000  
 Median :  2.00   Median : 0.0000   Median : 1.830   Median :1.000  
 Mean   : 11.99   Mean   : 0.1001   Mean   : 1.469   Mean   :1.218  
 3rd Qu.:  9.00   3rd Qu.: 0.0000   3rd Qu.: 2.440   3rd Qu.:2.000  
 Max.   :504.00   Max.   :25.0000   Max.   :11.000   Max.   :7.000

plot(census, "P25_TO_34", type="interval", breaks=5, plg=list(x="topright"))

Let’s do some interactive visualization with mapview.

library(mapview)
mapview(rstack$elevation) + census

Number of pixels is above 5e+05.Only about 5e+05 pixels will be shown.
You can increase the value of `maxpixels` to 2025000 to avoid this.

We highly recommend you to check the tmap package to make really appealing and publication ready maps.

To exemplify the use of write_* functions, let’s do a simple operation with the landuse raster map. We will apply a custom function that makes NULL all values less than 4.

result <- app(rstack$landuse, fun=function(x){ x[x < 4] <- NA; return(x)} )
plot(result)

write_RAST(result, "result_from_R", overwrite = TRUE)

Warning in execGRASS("r.in.gdal", flags = flags, input = tf, output = vname): The command:
r.in.gdal --overwrite input=/tmp/Rtmp6pJLgV/file218d3509b7c2a.grd output=result_from_R
produced at least one warning during execution:
WARNING: Raster map <result_from_R> already exists and will be overwritten
Importing raster map <result_from_R>...
   0%   3%   6%   9%  12%  15%  18%  21%  24%  27%  30%  33%  36%  39%  42%  45%  48%  51%  54%  57%  60%  63%  66%  69%  72%  75%  78%  81%  84%  87%  90%  93%  96%  99% 100%

WARNING: Raster map <result_from_R> already exists and will be overwritten
Importing raster map <result_from_R>...
   0%   3%   6%   9%  12%  15%  18%  21%  24%  27%  30%  33%  36%  39%  42%  45%  48%  51%  54%  57%  60%  63%  66%  69%  72%  75%  78%  81%  84%  87%  90%  93%  96%  99% 100%

SpatRaster read into GRASS using r.in.gdal from memory

execGRASS("g.list", parameters = list(type="raster", pattern="result*"))

result_from_R

B. You don’t have nor need GRASS projects

In the case you don’t have, need or care about GRASS projects, the initGRASS() function allows you to create temporary projects to use GRASS modules on R objects. This is equivalent to what QGIS does when you use GRASS tools via the Processing Toolbox.

So, the workflow is pretty similar to case A, except that we’ll create a temporary project based on the extent and CRS of a raster or vector R object. Hence, we do not pass project or mapset names but we need to pass a reference spatial grid. Then, we’ll write our R objects into the temporary GRASS project, run the desired processes, export the outputs back to R environment and done.

Let’s start with getting some spatial data, eg. a raster file, into R.

library(terra)

f <- system.file("ex/elev.tif", package="terra")
r <- rast(f)
plot(r)

Now, we’ll load the rgrass library and init GRASS with SG = r, so a GRASS project is internally created with r’s CRS (BTW, you can check that with crs(r)).

library(rgrass)

# path to GRASS binaries
grassbin <- system("grass --config path", intern = TRUE)

# start GRASS GIS from R
initGRASS(gisBase = grassbin, 
          home = tempdir(),
          SG = r, 
          override = TRUE)

gisdbase    /tmp/Rtmp6pJLgV 
location    file218d35570e7e6 
mapset      file218d33e5c586f 
rows        90 
columns     95 
north       50.19167 
south       49.44167 
west        5.741667 
east        6.533333 
nsres       0.008333333 
ewres       0.008333326 
projection:
 GEOGCRS["WGS 84",
    ENSEMBLE["World Geodetic System 1984 ensemble",
        MEMBER["World Geodetic System 1984 (Transit)"],
        MEMBER["World Geodetic System 1984 (G730)"],
        MEMBER["World Geodetic System 1984 (G873)"],
        MEMBER["World Geodetic System 1984 (G1150)"],
        MEMBER["World Geodetic System 1984 (G1674)"],
        MEMBER["World Geodetic System 1984 (G1762)"],
        MEMBER["World Geodetic System 1984 (G2139)"],
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ENSEMBLEACCURACY[2.0]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]

Now, we can write our SpatRaster into the GRASS temporary project.

write_RAST(r, "terra_elev")

Importing raster map <terra_elev>...
   0%   3%   6%  10%  13%  16%  20%  23%  26%  30%  33%  36%  40%  43%  46%  50%  53%  56%  60%  63%  66%  70%  73%  76%  80%  83%  86%  90%  93%  96% 100%

SpatRaster read into GRASS using r.in.gdal from file

Let’s check it is indeed within the project and run the GRASS module r.slope.aspect on it.

execGRASS("g.list", type = "raster")

terra_elev

execGRASS("r.slope.aspect", 
          elevation = "terra_elev", 
          slope = "slope",
          aspect = "aspect")

   0%   3%   6%   9%  12%  15%  18%  21%  24%  27%  30%  34%  37%  40%  43%  46%  49%  52%  56%  59%  62%  65%  68%  71%  74%  78%  81%  84%  87%  90%  93%  96% 100%
Aspect raster map <aspect> complete
Slope raster map <slope> complete

execGRASS("g.list", type = "raster")

aspect
slope
terra_elev

Let’s get slope and aspect maps into R

grass_maps <- read_RAST(c("aspect", "slope"))

Checking GDAL data type and nodata value...
   2%   5%   8%  11%  14%  17%  20%  23%  26%  30%  33%  36%  40%  43%  46%  50%  53%  56%  60%  63%  66%  70%  73%  76%  80%  83%  86%  90%  93%  96% 100%
Using GDAL data type <Float32>
Exporting raster data to RRASTER format...
   2%   5%   8%  11%  14%  17%  20%  23%  26%  30%  33%  36%  40%  43%  46%  50%  53%  56%  60%  63%  66%  70%  73%  76%  80%  83%  86%  90%  93%  96% 100%
r.out.gdal complete. File </tmp/Rtmp6pJLgV/file218d362249918.grd> created.
Checking GDAL data type and nodata value...
   2%   5%   8%  11%  14%  17%  20%  23%  26%  30%  33%  36%  40%  43%  46%  50%  53%  56%  60%  63%  66%  70%  73%  76%  80%  83%  86%  90%  93%  96% 100%
Using GDAL data type <Float32>
Exporting raster data to RRASTER format...
   2%   5%   8%  11%  14%  17%  20%  23%  26%  30%  33%  36%  40%  43%  46%  50%  53%  56%  60%  63%  66%  70%  73%  76%  80%  83%  86%  90%  93%  96% 100%
r.out.gdal complete. File </tmp/Rtmp6pJLgV/file218d36a6db337.grd> created.

grass_maps

class       : SpatRaster 
dimensions  : 90, 95, 2  (nrow, ncol, nlyr)
resolution  : 0.008333326, 0.008333333  (x, y)
extent      : 5.741667, 6.533333, 49.44167, 50.19167  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
sources     : file218d362249918.grd  
              file218d36a6db337.grd  
names       :       aspect,      slope 
min values  :   0.08974174, 0.01416342 
max values  : 360.00000000, 7.22943783

Now that the output maps are back into our R environment, we can plot them, do further analysis or write them into other raster formats.

plot(grass_maps)

writeRaster(grass_maps, "grass_maps.tif", overwrite=TRUE)

There’s yet another way in which you can use GRASS and R together, and it involves calling R from the GRASS terminal. In this way, rgrass will read all the environmental variables of the GRASS session, and you won’t need to use initGRASS. This would be similar to calling jupyter notebook within a GRASS active session. It goes more or less like this:

Open GRASS GIS
Type R or rstudio & in the GRASS GIS terminal
Load rgrass library with library(rgrass)
Use read_VECT(), read_RAST() to read data from GRASS into R
Access GRASS GIS modules and data through execGRASS()
Write data (back) to GRASS database with write_VECT() and write_RAST()

Starting GRASS GIS...

          __________  ___   __________    _______________
         / ____/ __ \/   | / ___/ ___/   / ____/  _/ ___/
        / / __/ /_/ / /| | \__ \\_  \   / / __ / / \__ \
       / /_/ / _, _/ ___ |___/ /__/ /  / /_/ // / ___/ /
       \____/_/ |_/_/  |_/____/____/   \____/___//____/

Welcome to GRASS GIS 8.3.0
GRASS GIS homepage:                      https://grass.osgeo.org
This version running through:            Bash Shell (/bin/bash)
Help is available with the command:      g.manual -i
See the licence terms with:              g.version -c
See citation options with:               g.version -x
If required, restart the GUI with:       g.gui wxpython
When ready to quit enter:                exit

Launching <wxpython> GUI in the background, please wait...
[Raster MASK present]
GRASS nc_basic_spm_grass7/PERMANENT:~ > R

R version 4.3.1 (2023-06-16) -- "Beagle Scouts"
Copyright (C) 2023 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(rgrass)
GRASS GIS interface loaded with GRASS version: GRASS 8.3.0 (2023)
and location: nc_basic_spm_grass7
>

Main funtions

Basic Usage: Choose your own adventure

A. You have GRASS projects

B. You don’t have nor need GRASS projects

References